// SBEB核心算法实现

/**
 * 将边进行聚类（分组），用于后续骨架抽取等图处理
 * 
 * @param {Array<[string, string]>} edges - 原始边集合，形式为 [fromId, toId]
 * @param {Array<{ id: string, lon: number, lat: number }>} nodes - 所有节点，包含地理坐标
 * @param {number} k - 聚类数目（默认 3）
 * @returns {Array<Array<[string, string]>>} clusteredEdges - 聚类后的边，每个子数组是一类边
 */
export function clusterEdges(edges, nodes, k, maxIterations) {
    // 为每条边构建特征向量：包含方向向量（单位长度） + 中点坐标
    const features = edges.map(([from, to]) => {
        const p1 = nodes.find(n => n.id == from);
        const p2 = nodes.find(n => n.id == to);

        // 计算方向向量并归一化
        const dx = p2.lon - p1.lon;
        const dy = p2.lat - p1.lat;
        const len = Math.sqrt(dx * dx + dy * dy) || 1e-6;
        const dir = [dx / len, dy / len];

        // 计算边的中点（lon, lat）
        const mid = [(p1.lon + p2.lon) / 2, (p1.lat + p2.lat) / 2];

        // 返回完整特征向量：[方向 x, 方向 y, 中点 lon, 中点 lat]
        return [...dir, ...mid];
    });

    // 对特征向量使用 KMeans 聚类（不依赖外部库）
    const assignments = kMeans(features, k, maxIterations);

    // 初始化聚类结果结构：k 个数组，每个数组对应一个聚类
    const clusteredEdges = Array.from({ length: k }, () => []);

    // 根据聚类分配结果，将边归入对应类别
    assignments.forEach((clusterIndex, i) => {
        clusteredEdges[clusterIndex].push(edges[i]);
    });

    return clusteredEdges;
}


/**
 * 使用 K-Means 聚类算法将数据划分为 k 个簇
 *
 * @param {Array<Array<number>>} data - 输入数据，每个元素是一个数值向量
 * @param {number} k - 聚类数目
 * @param {number} maxIterations - 最大迭代次数（默认 100）
 * @returns {Array<number>} assignments - 每个数据点所属的聚类编号
 */
function kMeans(data, k, maxIterations) {
    const dim = data[0].length; // 数据维度（如二维空间中的点 dim = 2）

    // 随机初始化 k 个质心
    const centroids = initializeCentroids(data, k);

    // assignments[i] 表示 data[i] 属于第几个簇（质心索引）
    const assignments = new Array(data.length).fill(0);

    for (let iter = 0; iter < maxIterations; iter++) {
        let changed = false;

        // 1. 将每个数据点分配到最近的质心
        for (let i = 0; i < data.length; i++) {
            const distances = centroids.map(c => euclideanDistance(data[i], c));
            const minIndex = distances.indexOf(Math.min(...distances));
            // 如果簇分配发生变化，则标记 changed
            if (assignments[i] !== minIndex) {
                assignments[i] = minIndex;
                changed = true;
            }
        }

        // 如果本轮没有任何分配变化，说明已经收敛，提前终止
        if (!changed) break; 

        // 2. 重新计算每个簇的质心（均值）
        const sums = Array.from({ length: k }, () => new Array(dim).fill(0));
        const counts = new Array(k).fill(0);  // 每个簇的点数

        data.forEach((vec, i) => {
            const cluster = assignments[i];
            counts[cluster]++;
            for (let d = 0; d < dim; d++) {
                sums[cluster][d] += vec[d];
            }
        });

        for (let j = 0; j < k; j++) {
            if (counts[j] === 0) continue; // 避免除以 0（可能出现空簇）
            for (let d = 0; d < dim; d++) {
                centroids[j][d] = sums[j][d] / counts[j]; // 更新为均值
            }
        }
    }

    return assignments;
}

function euclideanDistance(a, b) {
    return Math.sqrt(a.reduce((sum, val, i) => sum + (val - b[i]) ** 2, 0));
}

function initializeCentroids(data, k) {
    const shuffled = [...data].sort(() => 0.5 - Math.random());
    return shuffled.slice(0, k).map(vec => [...vec]); // 深拷贝
}

// 计算欧几里得距离
function distance(a, b) {
    const dx = a[0] - b[0];
    const dy = a[1] - b[1];
    return Math.sqrt(dx * dx + dy * dy);
}

/**
 * 使用 Kruskal 算法从一组点中构造最小生成树（MST），
 * 返回由最短边组成的连接所有点的无环连通图。
 *
 * @param {Array} points - 点集合，每个点为 [x, y] 数组。
 * @returns {Array} mst - 骨架线段集合，格式为 [[p1, p2], [p3, p4], ...]
 */
function kruskal(points) {
    // 并查集初始化，每个点自成一个集合
    const parent = Array.from({ length: points.length }, (_, i) => i);

    // 并查集查找函数，带路径压缩
    function find(i) {
        if (parent[i] !== i) parent[i] = find(parent[i]);
        return parent[i];
    }

    // 并查集合并函数，将 i 和 j 所在集合合并
    function union(i, j) {
        const ri = find(i);
        const rj = find(j);
        if (ri !== rj) parent[ri] = rj;
    }

    // 1. 构造所有可能点对之间的边及其距离
    const edges = [];
    for (let i = 0; i < points.length; i++) {
        for (let j = i + 1; j < points.length; j++) {
            edges.push([i, j, distance(points[i], points[j])]);
        }
    }
    // 2. 按距离从小到大排序边
    edges.sort((a, b) => a[2] - b[2]);

    // 3. Kruskal 主体：从最短边开始依次加入，不构成环即可加入 MST
    const mst = [];
    for (const [i, j, d] of edges) {
        if (find(i) !== find(j)) {
            union(i, j);
            mst.push([points[i], points[j]]); // 用线段表示骨架边
        }
    }
    return mst;
}


/**
 * 对骨架线段进行自适应加密（插值细分），使每段长度不超过 maxSegmentLength。
 *
 * @param {Array} edges - 原始骨架线段数组，每个元素为 [[x1, y1], [x2, y2]]
 * @param {number} maxSegmentLength - 最大允许的子段长度，默认值为 10
 * @returns {Array} densifiedEdges - 加密后的线段数组
 */
function densifySkeletonEdgesAdaptive(edges, maxSegmentLength) {
    const densifiedEdges = [];
    for (const [p1, p2] of edges) {
        // 计算当前线段长度
        const len = distance(p1, p2);
        // 计算需要切分成多少段（向上取整以保证每段不超过最大长度）
        const segments = Math.ceil(len / maxSegmentLength);
        // 使用线性插值将线段分成若干小段
        for (let i = 0; i < segments; i++) {
            const t1 = i / segments;
            const t2 = (i + 1) / segments;
            // 线性插值计算子段起点 (x1, y1) 和终点 (x2, y2)
            const x1 = p1[0] * (1 - t1) + p2[0] * t1;
            const y1 = p1[1] * (1 - t1) + p2[1] * t1;
            const x2 = p1[0] * (1 - t2) + p2[0] * t2;
            const y2 = p1[1] * (1 - t2) + p2[1] * t2;
            // 存入加密后的子线段
            densifiedEdges.push([[x1, y1], [x2, y2]]);
        }
    }
    return densifiedEdges;
}

/**
 * 从聚类后的边中提取骨架线段，用于边绑定与变形阶段。
 * 每个聚类内的边会被表示为中点集合，然后构造最小生成树（MST），再进行骨架加密。
 *
 * @param {Map} clusteredEdges - Map<clusterId, edges>，每个聚类包含若干边 [fromId, toId]
 * @param {Array} nodes - 节点数组，每个节点为 { id, lon, lat }
 * @returns {Object} skeletons - 每个聚类的骨架线段集合，格式为：
 *   { clusterId1: [[p1, p2], [p3, p4], ...], clusterId2: [...], ... }
 */
export function extractSkeleton(clusteredEdges, nodes, maxSegmentLength) {
    // 构建节点查找表
    const nodeMap = new Map(nodes.map(n => [n.id, n]));
    const skeletons = {};

    // 遍历每个聚类
    for (const [clusterId, edges] of clusteredEdges.entries()) {
         // 1. 将聚类中的每条边转换为几何中点，形成点集
        const midpoints = edges.map(([from, to]) => {
            const p1 = nodeMap.get(from);
            const p2 = nodeMap.get(to);
            return [(p1.lon + p2.lon) / 2, (p1.lat + p2.lat) / 2];
        });

        // 2. 使用 Kruskal 算法从中点集合中构造最小生成树（骨架的初步形式）
        const mst = kruskal(midpoints);
        // 3. 对骨架线段进行加密（细分），增强后续吸附效果
        const skeleton = densifySkeletonEdgesAdaptive(mst, maxSegmentLength); 
        // 4. 存储该聚类的骨架结构
        skeletons[clusterId] = skeleton;
    }

    return skeletons;
}


/**
 * 将每条边绑定到其对应聚类的骨架上，用于 Skeleton-Based Edge Bundling。
 * 绑定方式采用 KNN（k 近邻）和高斯加权，获得一个加权平均的骨架点作为“吸附目标”。
 *
 * @param {Array} nodes - 节点数组，每个节点为 { id, lon, lat }
 * @param {Object} skeletons - 每个聚类对应的骨架线段集合，结构为：
 *   { clusterId1: [[p1, p2], [p3, p4], ...], clusterId2: [...], ... }
 * @param {Map} clusteredEdges - 映射每个聚类 ID 到其边集合的 Map，
 *   每条边为 [fromId, toId]
 * @param {number} k - KNN 的 k 值，表示选取多少个最近骨架点参与加权，默认 8
 * @param {number} sigma - 高斯核函数的标准差，控制加权距离衰减，默认 0.1
 * @returns {Array} bindings - 每条边的绑定信息数组，每个元素为：
 *   { edge: [fromId, toId], midpoint: [x, y], nearest: [x, y], clusterId }
 */
export function computeBindings(nodes, skeletons, clusteredEdges, k = 8, sigma = 0.1) {
    // 构建节点映射表，便于根据 ID 查找节点
    const nodeMap = new Map(nodes.map(n => [n.id, n]));
    const bindings = [];

    // 1. 为每个聚类，构建骨架点集（去重）
    const skeletonPointsByCluster = {};
    for (const [clusterId, segments] of Object.entries(skeletons)) {
        const points = new Set();
        for (const [p1, p2] of segments) {
            points.add(p1.join(','));
            points.add(p2.join(','));
        }
        // 转回数字坐标数组 [[x, y], ...]
        skeletonPointsByCluster[clusterId] = [...points].map(p => p.split(',').map(Number));
    }

    /**
     * 给定一个边的中点，从骨架点中找到 k 个最近点，并使用高斯加权平均，得到最终绑定点
     *
     * @param {Array} midpoint - 当前边的中点 [x, y]
     * @param {Array} skeletonPoints - 当前聚类的骨架点集合 [[x, y], ...]
     * @returns {Array} nearestPoint - 加权后的绑定点 [x, y]
     */
    function findWeightedSkeletonMidpoint(midpoint, skeletonPoints) {
        if (!skeletonPoints || skeletonPoints.length === 0) {
            return midpoint;
        }
        // 计算所有骨架点与中点的欧式距离平方
        const distances = skeletonPoints.map(p => ({
            point: p,
            dist: (p[0] - midpoint[0]) ** 2 + (p[1] - midpoint[1]) ** 2
        }));
        // 按距离升序排序，取前 k 个最近点
        distances.sort((a, b) => a.dist - b.dist);
        const nearest = distances.slice(0, k);

        // 高斯加权平均坐标
        let sumX = 0, sumY = 0, totalWeight = 0;
        for (const { point, dist } of nearest) {
            const weight = Math.exp(-dist / (2 * sigma * sigma));
            sumX += point[0] * weight;
            sumY += point[1] * weight;
            totalWeight += weight;
        }

        return [sumX / totalWeight, sumY / totalWeight];
    }

    // 2. 遍历每个聚类中的边，计算绑定信息
    for (const [clusterId, clusterEdges] of clusteredEdges.entries()) {
        const skeletonPoints = skeletonPointsByCluster[clusterId];
        for (const [from, to] of clusterEdges) {
            // 获取边两端节点位置
            const p1 = nodeMap.get(from);
            const p2 = nodeMap.get(to);
            // 计算边的几何中点
            const midpoint = [(p1.lon + p2.lon) / 2, (p1.lat + p2.lat) / 2];
            // 找到加权后的绑定点
            const nearestPoint = findWeightedSkeletonMidpoint(midpoint, skeletonPoints);
            // 记录绑定信息
            bindings.push({
                edge: [from, to],
                midpoint,
                nearest: nearestPoint,
                clusterId
            });
        }
    }

    return bindings;
}


/**
 * 生成变形后的边路径（用于绘制 Polyline）
 * 利用骨架点作为“控制点”，通过二次 B 样条插值生成一条平滑曲线，
 * 实现 edge bundling 的视觉效果。
 *
 * @param {Array} bindings - 每条边的绑定信息数组，每个元素包含：
 *   - edge: [fromId, toId] 原始边
 *   - midpoint: 原始边的中点坐标 [lon, lat]
 *   - nearest: 该边绑定到的骨架点坐标 [lon, lat]
 *   - clusterId: 聚类 ID
 * @param {Array} nodes - 所有原始节点，结构为 { id, lon, lat }
 * @param {number} segments - 插值段数，决定曲线的平滑度，默认为 20
 * @param {number} alpha - 控制中点与骨架点之间的插值强度，范围 0~1，越大越靠近骨架点，默认为 0.9
 * @returns {Array} warpedEdges - 每条边变形后的曲线坐标数组，形如 [[p0, p1, ..., pn], ...]
 */
export function warpEdges(bindings, nodes, segments = 20, alpha = 0.9) {
    const nodeMap = new Map(nodes.map(n => [n.id, n]));
    const warpedEdges = [];

    for (const { edge: [from, to], nearest } of bindings) {
        const source = nodeMap.get(from);
        const target = nodeMap.get(to);
        const p0 = [source.lon, source.lat];
        const p2 = [target.lon, target.lat];

        const midpoint = [(p0[0] + p2[0]) / 2, (p0[1] + p2[1]) / 2];
        const control = [
            midpoint[0] * (1 - alpha) + nearest[0] * alpha,
            midpoint[1] * (1 - alpha) + nearest[1] * alpha
        ];

        const curve = generateBSpline(p0, control, p2, segments);
        warpedEdges.push(curve);
    }

    return warpedEdges;
}


/**
 * 使用二阶B样条对3个点生成一段光滑曲线
 * @param {Array} p0 - 起点
 * @param {Array} p1 - 控制点
 * @param {Array} p2 - 终点
 * @param {Number} segments - 插值段数
 * @returns {Array} 插值点数组
 */
function generateBSpline(p0, p1, p2, segments) {
    const result = [];
    for (let t = 0; t <= 1; t += 1 / segments) {
        const b0 = (1 - t) * (1 - t);
        const b1 = 2 * (1 - t) * t;
        const b2 = t * t;
        const x = b0 * p0[0] + b1 * p1[0] + b2 * p2[0];
        const y = b0 * p0[1] + b1 * p1[1] + b2 * p2[1];
        result.push([x, y]);
    }
    return result;
}
